Dynamical scaling of fragment distribution in drying paste 
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We reproduce patterns of drying paste by means of smoothed particle hydrodynamics which is the 
one of methods for solving the equations of continuum in the Lagrangian description. In addition to 
reproduce a realistic pattern, we find that average size of fragments decays in proportion to inverse 
time in the case of a linear drying process. Distributions of the size of the fragments are obtained 
depending on the time. We find a universal scaling distribution by scaling analysis with the average 
size of the fragment. 

PACS numbers: 46.50.+a, 62.20.mt, 89.75.Kd 



I. INTRODUCTION 

We often see patterns of fracture on drying lakes, 
paddy fields and so on. It is well-known that fractures 
of drying paste show different properties depending on 
the thickness of the paste. The drying cracks are clas- 
sified into two types. In the case that the thickness is 
larger than a horizontal length of system, the fragment's 
characteristic size on upper surface is limited and the 
cracks run slowly along depth direction. In this case, the 
cracks make prismatic structures. On the other hand, in 
the case that thickness is shorter than the size of system, 
the cracks approach the bottom of container immediately 
and do not make columnar structures. In this case, some 
interesting properties have been reported. Groismanp] 
did some experiments by using coffee granular and re- 
ported that the average size of fragments is proportional 
to the thickness of the paste. Furthermore, Nakahara 
and Matsuo[2|, [3| reported that some kind of paste re- 
members the force they received or the directions they 
flowed before drying on their crack patterns. Thus the 
properties of drying paste have been investigated vigor- 
ously, however, most of them were investigated only at 
the time when drying has finished. The patterns change 
every moment until the end of drying and the statistical 
properties change every moment. It is worth investigat- 
ing the dynamical and statistical properties for under- 
standing fractures of drying paste. There are studies for 
velocity of crack tips[4[ as the time-dependent properties, 
however, there are few studies of how patterns change 
with time. In this article, we focus on the properties of 
size of fragment as the time-dependent properties. 

Smoothed Particle Hydrodynamics (SPH) Q is the 
one of methods to calculate equations of continuum. 
SPH has been developed in the field of astrophysics to 
solve problems of compressible flow@, Q. Currently, 
SPH has been applied to calculations of incompressible 
fluids [9l [Toll and, furthermore, elastic [llj or visco-elastic 
materials |12]. There are also some studies of plastic- 
elastic materials [l3l [l4| and the formalization of treat- 
ments for ductile materials is nearly complete, however, 
the formalization of treatments for brittle materials have 
not been enough yet. 
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FIG. 1. Geometry of the paste we simulate in this paper. 
Thickness is H and length of a side is L. 



SPH belongs to a method in the Lagrangian descrip- 
tion. In the Lagrangian description, we do not have to 
make meshes. Therefore it is easy to calculate the equa- 
tions of systems which have complex free surfaces. In 
the case of drying fracture process, new free surfaces are 
created when cracks run. As a consequence, complex 
free surfaces are created. Therefore, SPH is a suitable 
method for simulations of drying paste. However, there 
is a problem in the implementation of brittle fractures in 
the SPH algorithm. Resolving this problem is the key for 
reproducing fractures of drying paste by using SPH. 



The aim of this article are (i) to reproduce a crack pat- 
tern of a thin drying paste with a continuum model by 
giving the implementation of the brittle material in the 
SPH algorithm, and (ii) to investigate the dynamical and 
statistical properties of size of fragment with the SPH 
simulation. To simplify the calculation, we make a two- 
dimensional continuum model from a three-dimensional 
Voigt visco-elastic continuum model. We investigate two 
characteristic quantities which are convenient for under- 
standing time evolution of the fragments size: the average 
size and the size distribution. As a result, we find that 
the average size decays in proportion to inverse time in 
the case of a linear drying process. In addition, we find 
a universal time-independent scaling distribution by the 
scaling analysis with the average size of the fragments. 
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II. MODEL OF DRYING PASTE 

In this paper, we take a thin layer of paste whose thick- 
ness is H and length of a side is L (See Fig.[T]). Let t?, ix, 
(j and p denote velocity, displacement, stress and density, 
respectively. In general, the equations of motion, conti- 
nuity and displacement of three dimensional continuum 
are described as follows: 



dv „ 



dp 
~dt 



-pV • v, 



du 
~dt 



(1) 



(2) 



(3) 



Here, d/dt is a time differential operator in the La- 
grangian description. In this paper, we treat the paste 
as the Voigt visco-elastic material. In addition, we take 
the effect of drying into account. The equation of stress 
a is given by 



& el 



T dry 



(4) 



Here, cr V i S is a viscous part. By using a strain velocity e 
given by 



{(V(g)i;) + (V(g)i;) T }, 



cr V je is described as 



rye, 



(5) 



(6) 



where 77 is a viscosity coefficient. a e i is an elastic part. 
By using a rotation velocity ft given by 



^ = ^{(V^^)-(V®v) T }, 
the time evolution of <j e \ is described as 



d(Tel 

dt 



ATr (c) I + 2pe + a el • ft - ft • a eh (8) 



where A and p are Lame's elastic constants. J is a unit 
tensor. 

(J dry is a drying part. Drying makes inner stress in- 
creased with time evolution. In general, we do not know 
the functional form of a dry- The stress is considered as 
the monotone increasing function of time t in initial stage 
of drying process. Therefore, we take (Tdry as a linear 
function of t. Furthermore, we assume that the effect of 
drying only appears in diagonal part of (Tdry In a pre- 
vious work[15|, the drying part was treated similarly to 
the present form. Eventually, (Tdry is described by 



® dry — V T tI , 



(9) 



where v T is the drying speed. 

To simplify the calculation, we approximate the equa- 
tion of motion to the two-dimensionalized one. In or- 
der to approximate, we use the discretization used in the 



study by Otsuki[16j. Let lower index i denote that the 
quantity is along the i direction, V{ at (x, y, 0) and 
<Ji Z at (x, y, H) must be satisfied boundary conditions as 
follows: 



a iz (x,y,H) = 0, 
Ui (x,y,0) = 0. 



(10) 
(11) 



If the thickness of paste H is thin enough, we can ig- 
nore the motions along z direction and consider that the 
space differential of a certain quantity along z direction 
is approximated to the difference between the quantities 
at the top (z = H) and the bottom (z = 0) of the paste. 
For example, let denote the strain along z direction 
and it at (x, y : 0) is given by 



(x,y,0) 



1 Ui (x,y,H) - Ui (x,y,0) 

2 H 

1 Uj (x,y,H) 

2 H 



(12) 



€i Z (x, y, 0) is also descritized similarly. When we consider 
the motion of the upper surface of the paste, we can 
divide the right-hand side of Eq. (pp) into the following 
equations: 



dvi 



dciiz 

~dz~ 



for ij = x,y at (x,y,H). (13) 



The second term in the right-hand side of Eq. ([T3]) works 
as a resistance force for the motion along i direction. This 
resistance force can be discretized as the following form: 

d(7iz ( m - Giz ^Ei y ' H ^ ~ Giz y ' °) 



H 

Giz (x, y, 0) 
H 

e iz (x,y,0) ei z (x,y,0) 
Tt V- 



H 



H 



(14) 

Here, we suppose the damping force, the second term of 
the right-hand side of Eq. (fT4]) , can be ignored because 
this term is effective only just in the time when cracks 
break out and is very smaller than the first term of the 
right-hand side of Eq. ([T4]) if the system shrinks slowly. 
Eventually, the two dimensional equation of motion is 
given by 



dv _ p 



(15) 



where u and a are two dimensional quantities, respec- 
tively. We calculate the motion of the upper surface by 
using Eqs. (O-© and (fT5]) . 

When we simulate fractures of continuum, we have 
to set a yield criterion. There are many yield criteri- 
ons which depend on the kind of simulated material. In 
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this paper, we use a criterion [17] by using local averaged 
stress a = Tr (a) /2, which corresponds to the pressure. 
Local averaged stress criterion is described as follows: (1) 
Calculate a at all positions in the material. (2) If a is 
grater than a definite yield stress cry , we make the stress 
at the position into zero. 



III. SMOOTHED PARTICLE 
HYDRODYNAMICS AND SIMULATION 
CONDITION 

SPHjil, El is the method for solving the equations of 
continuum by using a movable mesh point "particle" . All 
particles have the physical quantities. Let upper indexes 
written with capital letters denote particle numbers. Let 
r 1 denote the position of particle /. In SPH formula, a 
general physical quantity f 1 is given by 



lows: 



■h) 



m 



(16) 



where m J and p J are the mass and density of particle J. 
The summation of J is calculated by using all particles in 
the system. The function W(x; h) is a kernel function and 
has a positive parameter h called an "effective length". 
W(x\ h) is required to be the Dirac delta function under 
the limit of h \ 0. The accuracy of the simulation is 
related with the choice of the kernel function W(x; h). In 
previous works [llj of SPH simulation, a spline function 
is often used. In this paper, we choose the fifth order 
spline function as the kernel function as follows: 



W(x;h) = 



63 



478tt/i 2 



(17) 



where W5 (x) is the fifth order spline function: 



W 5 (x) 



with 



' q s (x) - Qq 2 (x) + 15qi (x) (\x\ < ±) 

q 3 (x)-6q 2 (x) (§<M<§) 

q 3 (x) (§<|*|<1) 

(KM), 



q k (x) = (k-3\x\f (A = 1,2, 3). 

In SPH description, a gradient of the physical quantity 
V f 1 is given by 



rn 



V/ 1 = ^f J VW (|r J -r J \-h) ^ 



(18) 



By using Eq. ([18]) and some kind of differential trans- 
portation, V • tr, V v and V • v are discretized as fol- 



u 



a 1 a J 



, (19) 



9 ) V 0* J ) 2 . 

(V ^f = ^ ^j- VW IJ (v J - v 1 ) , (20) 
J P 

(V-v) 1 = Tr ((V®vf) , (21) 



where VW IJ denotes VW (|r J — r J \ ; h) . 

When we calculate the time evolutions of Eqs. (j2j), (j3j), 
(|8j) and (jl3J) in SPH description, we do not have to calcu- 
late advection term since SPH is based on the Lagrangian 
description. Instead of calculating the advection term, 
we have to calculate the time evolutions of positions of 
particles as the following: 



dr 1 
~dt 



v 1 . 



(22) 



To avoid numerical instability, we use a velocity aver- 
aging method. In SPH simulations, interaction between 
particles has no repulsive core. Therefore, in extreme 
conditions, there are some possibilities of that particles 
close each other. It causes numerical instability of the 
calculation. To avoid this instability, the velocity averag- 
ing method is proposed by MonaghanpH [l9| . According 
to the Monaghan's description, the velocity of particle / 
is averaged as the following: 

v 1 ^ v 1 + eJ2 K - v ') W (\ rI ~ rJ \'i h )i 

(23) 

where p IJ denotes \ (// + p J ) and e is a tuning param- 
eter. In this paper, we choose e = 0.5. 

In the following study, Eqs. (O-© and (JTSJ) dis- 
cretized by using Eqs. (fT9|) ~ ([22]) and ([23]) are taken to 
be the basic equations of the model. 

We have to apply the local average stress criterion 
to SPH. It is plausible that fractures of drying paste 
are treated as brittle fractures. In simple consideration, 
we only have to make the stress of particle into zero, 
when the stress becomes grater than the yield stress cry • 
However, we can not reproduce brittle fractures in this 
way because stress can be relaxed rapidly; The stress 
on the surfaces of brittle cracks is relaxed as long as 
both the surfaces are within the length h. In SPH de- 
scription, there are some different treatments of brittle 
fractures [13, [H|. ^ n this paper, we achieve the brittle 
fracture by removing particles which satisfy local average 
stress criterion instead of resetting the stress. Using this 
method, the distance between two created surfaces be- 
comes enough large which prevents the stress relaxation. 

In the actual simulation, we prepare a square paste 
whose length of a side is L. We impose a periodic bound- 
ary condition on the paste. The values of parameters and 



Parameter 


Symbol 


Value 


First Lame's constant 


A 


1.0 


Second Lame's constant 


fi 


0.1 


Yield stress 


cry 


5.0xl0" 3 


Viscosity 


7] 


1.0 


Thickness of paste 


H 


0.316 


Drying speed 


V T 


2.2xlCT 5 


Length of a side 


T 

L 


iU.U 


Effective length 


h 


0.2 


Number of SPH particles 


N 


4.0xl0 4 


Intial density 


Po 


1.0 


Intial stress 




-0.01-O.Olxar 



TABLE I. Parameters and initial quantities used in the sim- 
ulation. 

initial physical quantities are summarized in TABLE HI 
These values are non-dimensionalized by po> A and rj. 
The unit of time and space are given by rj/ A and rj/y/poX. 
Initial positions of particles are located randomly to avoid 
being created anisotropic patterns. Initial displacement 
and velocity are taken to be zero. Initial density po is 
set the value shown in TABLE HI which is also the unit. 
Initial stress <xo is randomly chosen from the uniform dis- 
tribution shown in TABLE HI The mass m of the particle 
is calculated by using the initial density po and the initial 
position r by using the consistent condition: 

Po = Yl mJw (\ rI ~r J \;h). (24) 
J 

Equation ([24]) is linear simultaneous equations and can 
be solved as a large sparse matrix problem. In order to 
solve Eq. (|24|) . we use the conjugate gradient method. 
We calculate the time evolutions of quantities by us- 
ing the fourth order Runge-Kutta (RK4) and Predictor- 
Corrector (PC4) methods. RK4 is used both in a few 
initial time steps and in several time steps after removals 
of particles and PC4 is used in otherwise time steps. This 
is because the physical quantities become discontinuous 
when removals of particles occur and the calculation may 
become unstable if we only use PC4. 

IV. RESULTS 

A. Time evolution of a crack pattern 

Figure [2] shows the snapshots of a simulation of dry- 
ing model. From the upper-left figure, time passes in 
alphabetical order. In this figure, we can see some nucle- 
ations of cracks before the cracks run. There are some 
cases where the cracks do not grow. We have not under- 
stood the reason why there are cases whether cracks run 
or not. Cracks run almost straight. When a crack runs 
against into another one, they cross each other normally. 
As a result, fragments become almost convex polygonal 
shape. In the simulation, crack tips' velocities are the 
fastest at the time when the cracks break out and decay 













(a) t =30.0 


(b) t =60.0 


(c) t =120.0 










m 


(d) t =180.0 


(e) t =240.0 


(f) t =270.0 



FIG. 2. Snapshots of time evolution of drying model simu- 
lated with parameters shown in TABLE[H From the upper-left 
figure, time passes in alphabetical order. Color shows magni- 
tude of average stress. Black region is the crack. 




Time /; 

FIG. 3. Time evolution of average size (S)t' Vertical axis in- 
dicates reciprocal of {S)t- Horizontal axis indicates the time. 



with approaching to another one. These behaviors are 
very similar to real patterns of drying paste. 

Color shows a magnitude of the average stress normal- 
ized by the yield stress. In this figure, we can see that 
the average stress is almost zero along the edges of each 
fragments and becomes to maximum on the position near 
the center of area. We can consider that this is because 
fragments are almost convex polygonal shape. Since a 
fragment shrinks isotropically, if the fragment is convex 
polygonal shape, stress concentration tends to occur near 
the center of area and rarely occurs along the edges. 



B. Time evolution of average size of fragments 

In order to investigate the time evolution of average 
size of fragments, we prepare a hundred samples with 
a different initial stress. To calculate each area of frag- 
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FIG. 4. Size distributions of fragments: Horizontal axis indi- 
cates size of fragment and vertical axis indicates the distribu- 
tion. The distributions from t = 90.0 to t — 300.0 are shown 
in this figure. 



merits, we binarize the snapshots as following way: We 
discretized the snapshot specially with a square mesh. 
On each mesh point cc, a binary value (j) (x) is calculated 

by 



4>{x) 



for (j> < ]r 
otherwise, 



■h) 



where </>o is a threshold value. If enough particles sur- 
round the mesh point a?, 4>{x) becomes 1. This means 
that the position is the inside of a fragment. Conversely, 
if few particles surround the mesh point a?, (j) (x) becomes 
0. This means that the position is the outside of a frag- 
ment and a part of cracks. We can calculate the area 
of fragment by clustering and summing up the value of 
cj). In this paper, we choose 0o = 0.8 and the interval of 
square mesh Ax = 0.01A 

Let (S)t denote the average size of fragments at time 
t. Figure [3] shows the time evolution of (S)t- As we 
see in this figure, the reciprocal value of (S) t evolves in 
proportion to time except for the initial stage. In other 
words, (S)t evolves in inverse proportion to time in long 
time. In the initial stage, the time which is earlier than 
t ~ 150, we find that there are some cases of that an 
initial fragment have not broken into two or more pieces 
yet and ones of that an initial one have already broken. 

This property can be explained by means of an dimen- 
sional analysis. We consider an over-damped equation 
of motion which is made by ignoring the inertial term of 
Eq. ([T5]) . The equation is given by 



dcri 



- JL V . 
dx n ~ H 2 l ' 



(25) 



Furthermore, the constitutive equation ignoring the vis- 
cous term of Eq. dU and using a general drying stress 



F (t) is given by 



duj 
dxj 



duj_ 

dxi 



(26) 



We investigate Eqs. ([25]) and ([26]) by means of dimen- 
sional analysis. Let U (£), 5 (t) and L (t) denote a charac- 
teristic displacement, stress and length of the fragment, 
respectively. Note that the space differential operator 
d/dx is replaced into 1/L(t). By using these charac- 
teristic quantities, Eqs. ([25]) and ([26]) are rewritten as 
follows: 



L(t) 



3(t) = (A + 2 M ) 



U(t) 
L{t) 



F(t). 



(27) 



(28) 



From these equations, we obtain the relationship between 
E(t), L(t) and F (t) as 



L 2 (t) = 



H 



V F(t)/E(t)-r 



(29) 



The value of H (t) is limited by the yield stress. Therefore 
the time evolution of H (t) is ignored from the time evo- 
lution of L 2 (t). Average size (S)t must be corresponding 
to L 2 (t) which means a characteristic area. Eventually, 
(S)t is proportional to 1/F (t). Therefore (S) t is propor- 
tional to 1/t because F (t) is proportional to t in our case. 
As a practical matter, (S)t must be described by super- 
position of various L 2 (t), however, it does not affect the 
time dependence of (S) t in long time. 



C. Time evolution of size distributions of fragments 

Next we investigate the time evolution of size distribu- 
tions of fragments. Let P (S; t) denote the size distribu- 
tion at time t. Figure [4] shows the results. As we can see 
in this figure, the peak of distributions goes to smaller 
side with time evolution. This shift is trivial behavior. 
In the case of the present model, this peak will go to zero 
in the long-time limit since the end of drying is not taken 
into the account in this model. 

As a property of distributions, we find that the distri- 
butions can be scaled by their averages (see Fig. [5] and 
the details in Fig. [6]). In other words, by using a dimen- 
sionless variable X defined as 



X = S/{S) t 



(30) 



raw distributions P (5; t) can be transformed into a time- 
independent distribution P (X) as the following: 



(S)tP 



s 

JS)t 



■;t 



P(X). 



(31) 



However, this property does not hold in initial stage. It 
is clearly observed in Fig. [5] before t = 180.0. After the 
initial stage, Eq. ([3T]) holds very well with time. 
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t=180.0 
t=210.0 
t=240.0 
t=270.0 
t=300.0 



Scaled area 



0.5 1 1.5 2 

Scaled area J£ 



2.5 



FIG. 5. Size distributions scaled by their averages: Horizontal 
axis indicates scaled size X. Vertical axis indicates the scaled 
distribution. The distributions from t = 90.0 to t = 300.0 are 
shown in this figure. 



FIG. 6. This figure shows the scaled distributions of FIG [5] 
from t 180.0 to t 300.0. 



The initial stage of time evolution of (S)t corresponds 
the scaling violation stage. By using Eq. (|29|h we esti- 
mate the time scale of the initial stage. This time scale 
is determined by the time when F (t) /E (t) is greater 
than 1 in Eq. (^J). If F (t) /3 (t) is greater than 1, av- 
erage size decays in proportion to H (t) /F (£). As it has 
been mentioned above, the region of 5 (t) is limited by 
yield stress cry. Furthermore, the lower limit can be es- 
timated as greater than zero because the paste has been 
shrunk and received tension effectively. As a rough esti- 
mate, supposing that H (t) equals ay, the time scale of 
the initial stage is evaluated to oy jv T ~ 227. The actual 
time scale might be smaller this value, because S (t) is 
a characteristic stress in the fragments and smaller than 



V. CONCLUSIONS AND DISCUSSIONS 

In this article, we modeled a thin drying paste and re- 
produced crack patterns by means of the smoothed par- 
ticle hydrodynamics. Furthermore, we found two prop- 
erties: The power-law decaying of average size and the 
scaling law of size distributions by their averages. Ac- 
cording to a dimensional analysis, we found the relation- 
ship between the average size and the drying stress. 

In our model, since the average size is proportional to 
inverse time, the number of fragments increases in pro- 
portion to time. This result shows that the drying frac- 
ture process is not a simple dividing system such as a cell 
division whose number of cells increases exponentially. 
For understanding this result, it is necessary to consider 
the dividing process that the dynamics of drying fracture 
process is considered. Most of typical distributions, such 
as normal and log- normal distributions, have two charac- 



teristic parameters, an average and a variance, and they 
can be scaled by using two parameters. However, the 
distribution in our model can be scaled by only one pa- 
rameter, the average. Therefore, we can consider that 
these distributions are unknown distributions which can 
not be described by some typical distributions. We have 
not been able to identify the functional form of these dis- 
tributions yet. Curz et al. [20] have shown that a mass 
distribution is scaled by their average mass in a simu- 
lation of fragmentation of hard-core granular gases with 
different restitution coefficients. It does not have direct 
relation with our results. But there are some analogous 
properties with the present scaled distribution. 

The patterns of our model are similar to the experi- 
mental patterns of drying paste. Therefore, we are re- 
quired to develop how to measure the similarity of pat- 
terns quantitatively and confirm the similarity with the 
experiments. We have confirmed only the power decay- 
ing and the scaling law of the distribution in experimental 
results preliminarily. [2lj In order to give a precise con- 
clusion of both properties, however, we must obtain and 
analyze a lot of experimental data. Size distributions of 
fragments can be scaled by their averages. So if we could 
know the behavior of average size of fragments, we can 
know the future distributions of size from a initial distri- 
bution. If actual experiments have this scaling law, it is 
possible to predict the distributions. 
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